

#include "slep.h"
#include "hr_time.h"
#include "epph.h"
#include "float.h"
#include "Logistic.h"

FUNVAL mtLogisticRm1(double* A, double* y, double z, OPTS opt);
FUNVAL mtLogisticRm2(double* A, double* y, double z, OPTS opt);
FUNVAL mtLogisticRm3(double* A, double* y, double z, OPTS opt);


FUNVAL mtLogisticR(double* A, double* y, double z, OPTS opt){
	if(opt.mFlag== 0 && opt.lFlag == 0 )
		return mtLogisticRm1(A,y,z,opt);
	else if(opt.mFlag == 1 && opt.lFlag == 0 && opt.q == 2)
		return mtLogisticRm2(A,y,z,opt);
	else if(opt.mFlag == 1 && opt.lFlag == 1 && opt.q == 2)
		return mtLogisticRm3(A,y,z,opt);
	else
		errorOut("Please set the values of opt.mFlag and opt.lFlag to correct values before executing mtLogisticR.");
}


FUNVAL mtLogisticRm1(double* A, double* y,double z,OPTS opt){
	int m,n,k, iterStep,i,j, *ind, bFlag;
	double q, lambda, q_bar, lambda_max;
	double *m1, *m2;
	int *p_flag;
	double *ATb, *x, *c, *Ax;
	double *cp, *ccp, *xp, *xxp, *Axp;
	double norm_xp, norm_xxp, *funVal, *ValueL;
	double fVal, fValp;
	double alpha, alphap, L, beta;
	double fun_x, fun_s, l_sum, r_sum;
	double *gc, *sc;
	double *g, *v, *aa, *bb, *s, *As, *tt, *b, *prob;
	FUNVAL ret;

	m = opt.m;
	n = opt.n;
	if(z<= 0)
		errorOut("mtLogisticR: z should be positive");
	// sll_opts

	if(!hasGot(&opt,"ind"))
		errorOut("mtLogisticR: ind should be specified");
	else{
		ind = opt.ind;
		k = opt.ind_size -1;
		if(ind[k] != m)
			errorOut("mtLogisticR: Check opt.ind");
	}
	if( !hasGot(&opt,"q"))
		q = opt.q = 2;
	else if( (q=opt.q) <1)
		errorOut("mtLogisticR: q should be greater than 1");

	initNormalization(&opt,A);
	m1 = (double *) malloc( sizeof(double) * k);
	m2 = (double *) malloc( sizeof(double) * k);
	p_flag = (int *) malloc( sizeof(int) * m);
	b = (double *) malloc( sizeof(double) * m);
	ATb = (double *) malloc( sizeof(double) * n*k);
	x = (double *) malloc(sizeof(double) * n*k);
	c = (double *) malloc(sizeof(double) * k);
	Ax = (double *) malloc(sizeof(double) * m);
	cp = (double *) malloc(sizeof(double) * k);
	ccp = (double *) malloc(sizeof(double) * k);
	gc = (double *) malloc(sizeof(double) * k);
	sc = (double *) malloc(sizeof(double) * k);
	xp = (double *) malloc(sizeof(double) * n*k);
	xxp = (double *) malloc(sizeof(double) * n*k);
	Axp = (double *) malloc(sizeof(double) * m);
	g = (double *) malloc(sizeof(double) * n*k);
	v = (double *) malloc(sizeof(double) * n*k);
	aa = (double *) malloc(sizeof(double) * m);
	bb = (double *) malloc(sizeof(double) * m); 
	prob = (double *) malloc(sizeof(double) * m); 
	s = (double *) malloc(sizeof(double) * n*k);
	As = (double *) malloc(sizeof(double) * m);
	tt = (double *) malloc(sizeof(double) * maxof2(n,m));
	ValueL = (double *) malloc(sizeof(double) * opt.maxIter);
	funVal = (double *) malloc(sizeof(double) * opt.maxIter);


	memset(m1,0,sizeof(double)*k);
	memset(m2,0,sizeof(double)*k);
	for(i=0,j=0;i<m;i++){
		p_flag[i] = ((y[i]>0.9999 && y[i]<1.0001)?1:0);
		if(i >= ind[j+1])
			j++;
		m1[j] += p_flag[i];
		m2[j] = (ind[j+1] - ind[j]) - m1[j];
	}

	if(opt.rFlag == 0)
		lambda = z;
	else if( z<0 || z>1)
		errorOut("mtLogisticR: opt.rFlag=1 and z should be [0,1]");
	else{
		for(i=0,j=0;i<m;i++){
			if(i >= ind[j+1])
				j++;
			b[i] = (p_flag[i] * m1[j]/(m1[j]+m2[j]))
				- ( (1^p_flag[i]) * m2[j]/(m1[j]+m2[j]));
		}
		
		// ATb = zeros(n,k);
		memset(ATb,0,sizeof(double)*n*k);

		// Compute ATb
		for(i=0;i<k;i++){
			normalizedmv_lowlevel('T',ind[i+1]-ind[i] ,n,&A[ind[i]],m,&b[ind[i]],tt,ind[i],0, &opt);
			dcopy(n,tt,1,&ATb[i*n],1);
		}

		if(q==1)
			q_bar = DBL_MAX;
		else if ( q> 1.0e6 )
			q_bar = 1.0;
		else
			q_bar = q/(q-1);
		lambda_max = 0.0;
		for(i=0; i<n ; i++)
			lambda_max = maxof2(lambda_max ,dnrm(k,&ATb[i],n,q_bar));
		lambda_max = lambda_max / m;
		lambda = z * lambda_max;
	}

	if( opt.init == 2){
		memset(x,0,sizeof(double)*n*k);
		memset(c,0,sizeof(double)*k);
	}
	else{
		if(hasGot(&opt,"x0"))
			x = opt.x0;
		else
			memset(x,0,sizeof(double)*n*k);
		if( hasGot(&opt,"mt_c0"))
			c = opt.mt_c0;
		else // c = log(m1./m2);
			for(i=0;i<k;i++)
				c[i] = log(m1[i]/m2[i]);
	}

	for(i=0; i < k; i++)
		normalizedmv_lowlevel('N',ind[i+1]-ind[i],n,&A[ind[i]],m,&x[i*n],&Ax[ind[i]],ind[i],0, &opt);

	if( opt.mFlag == 0 && opt.lFlag == 0){
		bFlag = 0;

		L = 1.0l/m;

		//xp=x; Axp=Ax; xxp=zeros(n,k);
		//cp=c; ccp=zeros(1,k);
		dcopy(n*k,x,1,xp,1);
		dcopy(m,Ax,1,Axp,1);
		memset(xxp,0,sizeof(double)*n*k);
		dcopy(k,c,1,cp,1);
		memset(ccp,0,sizeof(double)*k);

		alphap = 0; 
		alpha=1;

		for(iterStep=0;iterStep<opt.maxIter;iterStep++){
			beta = (alphap-1)/alpha;
			// s = x + beta* xxp;
			dcopy(n*k,x,1,s,1);
			daxpy(n*k,beta,xxp,1,s,1);

			//sc = c + beta * ccp;
			dcopy(k,c,1,sc,1);
			daxpy(k,beta,ccp,1,sc,1);

			// As = Ax + beta * (Ax - Axp);
			dcopy(m,Ax,1,As,1);
			dscal(m,1+beta,As,1);
			daxpy(m,-beta,Axp,1,As,1);

			dcopy(m,As,1,aa,1);
			//aa= - diag(y) * (A * s + sc)
			for(i=0; i < k; i++)
				daxpy(ind[i+1]-ind[i],sc[i],oneVector,1,&aa[ind[i]],1);
			dtbmv('U','N','N',m,0,y,1,aa,1);
			dscal(m,-1.0,aa,1);

			fun_s = 0;
			for(i=0; i<m;i++){
				bb[i] = maxof2(aa[i],0.0);
				fun_s += ( log(exp(-bb[i])+exp(aa[i]-bb[i])) + bb[i]);
			}
			fun_s = fun_s / m;

			//b = -y .* (1-prob);
			for(i=0;i<m;i++){
				prob[i] = 1.0l/(1.0l+exp(aa[i]));
				b[i] = (- y[i] * (1.0l - prob[i])) / m;
			}
			

//			memset(gc,0,sizeof(double)*k);
			for(i=0;i<k; i++)
				gc[i] = ddot(ind[i+1]-ind[i],&b[ind[i]],1,oneVector,1);

			// Compute ATb
			for(i=0;i<k;i++){
				normalizedmv_lowlevel('T',ind[i+1]-ind[i] ,n,&A[ind[i]],m,&b[ind[i]],tt,ind[i],0, &opt);
				dcopy(n,tt,1,&g[i*n],1);
			}

			//xp=x; Axp=Ax;
			//cp=c;
			dcopy(n*k,x,1,xp,1);
			dcopy(m,Ax,1,Axp,1);
			dcopy(k,c,1,cp,1);

			for(;;){
				//v = s - g/L; c = sc - gc/L;
				dcopy(n*k,s,1,v,1);
				daxpy(n*k,-1/L,g,1,v,1);
				dcopy(k,sc,1,c,1);
				daxpy(k,-1/L,gc,1,c,1);

				eppMatrix(x,v,n,k,lambda/L,q);

				//v = x-s;
				dcopy(n*k,x,1,v,1);
				daxpy(n*k,-1,s,1,v,1);

				// Compute Ax
				for(i=0; i < k; i++)
					normalizedmv_lowlevel('N',ind[i+1]-ind[i],n,&A[ind[i]],m,&x[i*n],&Ax[ind[i]],ind[i],0, &opt);

				dcopy(m,Ax,1,aa,1);
				//aa= - diag(y) * (A * x + c)
				for(i=0; i < k; i++)
					daxpy(ind[i+1]-ind[i],c[i],oneVector,1,&aa[ind[i]],1);
				dtbmv('U','N','N',m,0,y,1,aa,1);
				dscal(m,-1.0,aa,1);

				fun_x = 0;
				for(i=0; i<m;i++){
					bb[i] = maxof2(aa[i],0.0);
					fun_x += ( log(exp(-bb[i])+exp(aa[i]-bb[i])) + bb[i]);
				}
				fun_x = fun_x / m;

				r_sum = (ddot(n*k,v,1,v,1) + ddot(k,c,1,c,1) 
					+ ddot(k,sc,1,sc,1) - 2*ddot(k,c,1,sc,1))/2; 
				l_sum = fun_x - fun_s - ddot(n*k,v,1,g,1) 
					- ddot(k,c,1,gc,1) + ddot(k,sc,1,gc,1);

				if(r_sum <= 1e-20){
					bFlag=1;
					break;
				}

				if(l_sum <= r_sum * L)
					break;
				else
					L = maxof2( 2*L , l_sum/r_sum );
			}

			alphap = alpha;
			alpha = (1+ sqrt(4*alpha*alpha + 1))/2;

			ValueL[iterStep] = L;

			//xxp=x-xp;   ccp=c-cp;
			dcopy(n*k,x,1,xxp,1);
			daxpy(n*k,-1,xp,1,xxp,1);
			dcopy(k,c,1,ccp,1);
			daxpy(k,-1,cp,1,ccp,1);

			fVal = fun_x;
			for(i=0;i<n;i++)
				fVal += lambda* dnrm(k,&x[i],n,q);
			fValp = funVal[iterStep] = fVal;

			if(iterStep!=0)
				fValp = funVal[iterStep-1];

			if(bFlag)
				break;
			norm_xp = dnrm2(n*k,xp,1);
			norm_xxp = dnrm2(n*k,xxp,1);
			if(terminationCondition(&opt,fVal,fValp,x,norm_xxp,norm_xp,iterStep))
				break;
		}
	}
	ret.errorCode = 0;
	ret.funVal = funVal;
	ret.totIter = iterStep;
	ret.ValueL = ValueL;
	ret.x = x;
	return ret;
}

FUNVAL mtLogisticRm2(double* A, double* y,double z,OPTS opt){
	int m,n,k, iterStep,i,j, *ind, bFlag;
	double q, lambda, q_bar, lambda_max;
	double *m1, *m2;
	int *p_flag;
	double *ATb, *x, *c, *Ax;
	double *cp, *ccp, *xp, *xxp, *Axp;
	double norm_xp, norm_xxp, *funVal, *ValueL;
	double fVal, fValp;
	double alpha, alphap, L, beta;
	double fun_x, fun_s, l_sum, r_sum;
	double *gc, *sc;
	double *g, *v, *aa, *bb, *s, *As, *tt, *b, *prob;
	double *t, *tp;
	double *s_t, *u;
	FUNVAL ret;

	m = opt.m;
	n = opt.n;
	if(z<= 0)
		errorOut("mtLogisticR: z should be positive");
	// sll_opts

	if(!hasGot(&opt,"ind"))
		errorOut("mtLogisticR: ind should be specified");
	else{
		ind = opt.ind;
		k = opt.ind_size -1;
		if(ind[k] != m)
			errorOut("mtLogisticR: Check opt.ind");
	}
	if( !hasGot(&opt,"q"))
		q = opt.q = 2;
	else if( (q=opt.q) <1)
		errorOut("mtLogisticR: q should be greater than 1");

	initNormalization(&opt,A);
	m1 = (double *) malloc( sizeof(double) * k);
	m2 = (double *) malloc( sizeof(double) * k);
	p_flag = (int *) malloc( sizeof(int) * m);
	b = (double *) malloc( sizeof(double) * m);
	ATb = (double *) malloc( sizeof(double) * n*k);
	x = (double *) malloc(sizeof(double) * n*k);
	c = (double *) malloc(sizeof(double) * k);
	Ax = (double *) malloc(sizeof(double) * m);
	cp = (double *) malloc(sizeof(double) * k);
	ccp = (double *) malloc(sizeof(double) * k);
	gc = (double *) malloc(sizeof(double) * k);
	sc = (double *) malloc(sizeof(double) * k);
	xp = (double *) malloc(sizeof(double) * n*k);
	xxp = (double *) malloc(sizeof(double) * n*k);
	Axp = (double *) malloc(sizeof(double) * m);
	g = (double *) malloc(sizeof(double) * n*k);
	v = (double *) malloc(sizeof(double) * n*k);
	u = (double *) malloc(sizeof(double) * n*k);
	aa = (double *) malloc(sizeof(double) * m);
	bb = (double *) malloc(sizeof(double) * m); 
	prob = (double *) malloc(sizeof(double) * m); 
	s = (double *) malloc(sizeof(double) * n*k);
	As = (double *) malloc(sizeof(double) * m);
	tt = (double *) malloc(sizeof(double) * maxof2(n,m));
	ValueL = (double *) malloc(sizeof(double) * opt.maxIter);
	funVal = (double *) malloc(sizeof(double) * opt.maxIter);
	t = (double *) malloc(sizeof(double) * n);
	tp = (double *) malloc(sizeof(double) * n);
	s_t = (double *) malloc(sizeof(double) * n);


	memset(m1,0,sizeof(double)*k);
	memset(m2,0,sizeof(double)*k);
	for(i=0,j=0;i<m;i++){
		p_flag[i] = ((y[i]>0.9999 && y[i]<1.0001)?1:0);
		if(i >= ind[j+1])
			j++;
		m1[j] += p_flag[i];
		m2[j] = (ind[j+1] - ind[j]) - m1[j];
	}

	if(opt.rFlag == 0)
		lambda = z;
	else if( z<0 || z>1)
		errorOut("mtLogisticR: opt.rFlag=1 and z should be [0,1]");
	else{
		for(i=0,j=0;i<m;i++){
			if(i >= ind[j+1])
				j++;
			b[i] = (p_flag[i] * m1[j]/(m1[j]+m2[j]))
				- ( (1^p_flag[i]) * m2[j]/(m1[j]+m2[j]));
		}
		
		// ATb = zeros(n,k);
		memset(ATb,0,sizeof(double)*n*k);

		// Compute ATb
		for(i=0;i<k;i++){
			normalizedmv_lowlevel('T',ind[i+1]-ind[i] ,n,&A[ind[i]],m,&b[ind[i]],tt,ind[i],0, &opt);
			dcopy(n,tt,1,&ATb[i*n],1);
		}

		if(q==1)
			q_bar = DBL_MAX;
		else if ( q> 1.0e6 )
			q_bar = 1.0;
		else
			q_bar = q/(q-1);
		lambda_max = 0.0;
		for(i=0; i<n ; i++)
			lambda_max = maxof2(lambda_max ,dnrm(k,&ATb[i],n,q_bar));
		lambda_max = lambda_max / m;
		lambda = z * lambda_max;
	}

	if( opt.init == 2){
		memset(x,0,sizeof(double)*n*k);
		memset(c,0,sizeof(double)*k);
	}
	else{
		if(hasGot(&opt,"x0"))
			x = opt.x0;
		else
			memset(x,0,sizeof(double)*n*k);
		if( hasGot(&opt,"mt_c0"))
			c = opt.mt_c0;
		else // c = log(m1./m2);
			for(i=0;i<k;i++)
				c[i] = log(m1[i]/m2[i]);
	}

	for(i=0; i < k; i++)
		normalizedmv_lowlevel('N',ind[i+1]-ind[i],n,&A[ind[i]],m,&x[i*n],&Ax[ind[i]],ind[i],0, &opt);

	if( opt.mFlag == 1 && opt.lFlag == 0 && opt.q==2){
		bFlag = 0;

		L = 1.0l/m;

		//xp=x; Axp=Ax; xxp=zeros(n,k);
		//cp=c; ccp=zeros(1,k);
		dcopy(n*k,x,1,xp,1);
		dcopy(m,Ax,1,Axp,1);
		memset(xxp,0,sizeof(double)*n*k);
		dcopy(k,c,1,cp,1);
		memset(ccp,0,sizeof(double)*k);

		//t=sqrt(sum(x.^2, 2)); tp=t;
		for(i=0; i<n ; i++)
			t[i] = dnrm2(k,&x[i],n);
		dcopy(n,t,1,tp,1);


		alphap = 0; 
		alpha=1;

		for(iterStep=0;iterStep<opt.maxIter;iterStep++){
			beta = (alphap-1)/alpha;
			// s = x + beta* xxp;
			dcopy(n*k,x,1,s,1);
			daxpy(n*k,beta,xxp,1,s,1);

			//sc = c + beta * ccp;
			dcopy(k,c,1,sc,1);
			daxpy(k,beta,ccp,1,sc,1);

			// s_t= t + beta * (t-tp);
			dcopy(n,t,1,s_t,1);
			dscal(n,1.0 + beta, s_t,1);
			daxpy(n,-beta,tp,1,s_t,1);

			// As = Ax + beta * (Ax - Axp);
			dcopy(m,Ax,1,As,1);
			dscal(m,1+beta,As,1);
			daxpy(m,-beta,Axp,1,As,1);

			dcopy(m,As,1,aa,1);
			//aa= - diag(y) * (A * s + sc)
			for(i=0; i < k; i++)
				daxpy(ind[i+1]-ind[i],sc[i],oneVector,1,&aa[ind[i]],1);
			dtbmv('U','N','N',m,0,y,1,aa,1);
			dscal(m,-1.0,aa,1);

			fun_s = 0;
			for(i=0; i<m;i++){
				bb[i] = maxof2(aa[i],0.0);
				fun_s += ( log(exp(-bb[i])+exp(aa[i]-bb[i])) + bb[i]);
			}
			fun_s = fun_s / m;

			//b = -y .* (1-prob);
			for(i=0;i<m;i++){
				prob[i] = 1.0l/(1.0l+exp(aa[i]));
				b[i] = (- y[i] * (1.0l - prob[i])) / m;
			}
			

//			memset(gc,0,sizeof(double)*k);
			for(i=0;i<k; i++)
				gc[i] = ddot(ind[i+1]-ind[i],&b[ind[i]],1,oneVector,1);

			// Compute ATb
			for(i=0;i<k;i++){
				normalizedmv_lowlevel('T',ind[i+1]-ind[i] ,n,&A[ind[i]],m,&b[ind[i]],tt,ind[i],0, &opt);
				dcopy(n,tt,1,&g[i*n],1);
			}

			//xp=x; Axp=Ax;
			//cp=c;
			dcopy(n*k,x,1,xp,1);
			dcopy(m,Ax,1,Axp,1);
			dcopy(k,c,1,cp,1);

			//tp=t;
			dcopy(n,t,1,tp,1);

			for(;;){
				//u = s - g/L; 
				dcopy(n*k,s,1,u,1);
				daxpy(n*k,-1/L,g,1,u,1);

				//c = sc - gc/L;
				dcopy(k,sc,1,c,1);
				daxpy(k,-1/L,gc,1,c,1);

				//v= s_t - lambda / L;
				dcopy(n,s_t,1,v,1);
				daxpy(n,-lambda/L,oneVector,1,v,1);


				//eppMatrix(x,v,n,k,lambda/L,q);
				ep21R(x,t,u,v,n,k);

				//v = x-s;
				dcopy(n*k,x,1,v,1);
				daxpy(n*k,-1,s,1,v,1);

				// Compute Ax
				for(i=0; i < k; i++)
					normalizedmv_lowlevel('N',ind[i+1]-ind[i],n,&A[ind[i]],m,&x[i*n],&Ax[ind[i]],ind[i],0, &opt);

				dcopy(m,Ax,1,aa,1);
				//aa= - diag(y) * (A * x + c)
				for(i=0; i < k; i++)
					daxpy(ind[i+1]-ind[i],c[i],oneVector,1,&aa[ind[i]],1);
				dtbmv('U','N','N',m,0,y,1,aa,1);
				dscal(m,-1.0,aa,1);

				fun_x = 0;
				for(i=0; i<m;i++){
					bb[i] = maxof2(aa[i],0.0);
					fun_x += ( log(exp(-bb[i])+exp(aa[i]-bb[i])) + bb[i]);
				}
				fun_x = fun_x / m;

				r_sum = (ddot(n*k,v,1,v,1) + ddot(k,c,1,c,1) 
					+ ddot(k,sc,1,sc,1) - 2*ddot(k,c,1,sc,1) -2*ddot(n,t,1,s_t,1) + ddot(n,t,1,t,1) + ddot(n,s_t,1,s_t,1))/2 ; 
				l_sum = fun_x - fun_s - ddot(n*k,v,1,g,1) 
					- ddot(k,c,1,gc,1) + ddot(k,sc,1,gc,1);

				if(r_sum <= 1e-20){
					bFlag=1;
					break;
				}

				if(l_sum <= r_sum * L)
					break;
				else
					L = maxof2( 2*L , l_sum/r_sum );
			}

			alphap = alpha;
			alpha = (1+ sqrt(4*alpha*alpha + 1))/2;

			ValueL[iterStep] = L;

			//xxp=x-xp;   ccp=c-cp;
			dcopy(n*k,x,1,xxp,1);
			daxpy(n*k,-1,xp,1,xxp,1);
			dcopy(k,c,1,ccp,1);
			daxpy(k,-1,cp,1,ccp,1);

			fVal = fun_x + lambda*ddot(n,t,1,oneVector,1);

			fValp = funVal[iterStep] = fVal;

			if(iterStep!=0)
				fValp = funVal[iterStep-1];

			if(bFlag)
				break;
			norm_xp = dnrm2(n*k,xp,1);
			norm_xxp = dnrm2(n*k,xxp,1);
			if(terminationCondition(&opt,fVal,fValp,x,norm_xxp,norm_xp,iterStep))
				break;
		}
	}
	ret.errorCode = 0;
	ret.funVal = funVal;
	ret.totIter = iterStep;
	ret.ValueL = ValueL;
	ret.x = x;
	return ret;
}

FUNVAL mtLogisticRm3(double* A, double* y,double z,OPTS opt){
	int m,n,k, iterStep,i,j, *ind, bFlag;
	double q, lambda, q_bar, lambda_max;
	double *m1, *m2;
	int *p_flag;
	double *ATb, *x, *c, *Ax;
	double *cp, *ccp, *xp, *xxp, *Axp;
	double norm_xp, norm_xxp, *funVal, *ValueL;
	double fVal, fValp;
	double alpha, alphap, L, beta;
	double fun_x, fun_s, l_sum, r_sum;
	double *gc, *sc;
	double *g, *v, *aa, *bb, *s, *As, *tt, *b, *prob;
	double *t, *tp;
	double *s_t, *u;
	double gamma, *xnew, *cnew, *tnew, *Axnew, tao;
	FUNVAL ret;

	m = opt.m;
	n = opt.n;
	if(z<= 0)
		errorOut("mtLogisticR: z should be positive");
	// sll_opts

	if(!hasGot(&opt,"ind"))
		errorOut("mtLogisticR: ind should be specified");
	else{
		ind = opt.ind;
		k = opt.ind_size -1;
		if(ind[k] != m)
			errorOut("mtLogisticR: Check opt.ind");
	}
	if( !hasGot(&opt,"q"))
		q = opt.q = 2;
	else if( (q=opt.q) <1)
		errorOut("mtLogisticR: q should be greater than 1");

	initNormalization(&opt,A);
	m1 = (double *) malloc( sizeof(double) * k);
	m2 = (double *) malloc( sizeof(double) * k);
	p_flag = (int *) malloc( sizeof(int) * m);
	b = (double *) malloc( sizeof(double) * m);
	ATb = (double *) malloc( sizeof(double) * n*k);
	x = (double *) malloc(sizeof(double) * n*k);
	c = (double *) malloc(sizeof(double) * k);
	Ax = (double *) malloc(sizeof(double) * m);
	cp = (double *) malloc(sizeof(double) * k);
	ccp = (double *) malloc(sizeof(double) * k);
	gc = (double *) malloc(sizeof(double) * k);
	sc = (double *) malloc(sizeof(double) * k);
	xp = (double *) malloc(sizeof(double) * n*k);
	xxp = (double *) malloc(sizeof(double) * n*k);
	Axp = (double *) malloc(sizeof(double) * m);
	g = (double *) malloc(sizeof(double) * n*k);
	v = (double *) malloc(sizeof(double) * n*k);
	u = (double *) malloc(sizeof(double) * n*k);
	aa = (double *) malloc(sizeof(double) * m);
	bb = (double *) malloc(sizeof(double) * m); 
	prob = (double *) malloc(sizeof(double) * m); 
	s = (double *) malloc(sizeof(double) * n*k);
	As = (double *) malloc(sizeof(double) * m);
	tt = (double *) malloc(sizeof(double) * maxof2(n,m));
	ValueL = (double *) malloc(sizeof(double) * opt.maxIter);
	funVal = (double *) malloc(sizeof(double) * opt.maxIter);
	t = (double *) malloc(sizeof(double) * n);
	tp = (double *) malloc(sizeof(double) * n);
	s_t = (double *) malloc(sizeof(double) * n);
	xnew = (double *) malloc(sizeof(double) * n*k);
	tnew = (double *) malloc(sizeof(double) * n);
	cnew = (double *) malloc(sizeof(double) * k);
	Axnew = (double *) malloc(sizeof(double) * m);


	memset(m1,0,sizeof(double)*k);
	memset(m2,0,sizeof(double)*k);
	for(i=0,j=0;i<m;i++){
		p_flag[i] = ((y[i]>0.9999 && y[i]<1.0001)?1:0);
		if(i >= ind[j+1])
			j++;
		m1[j] += p_flag[i];
		m2[j] = (ind[j+1] - ind[j]) - m1[j];
	}

	if(opt.rFlag == 0)
		lambda = z;
	else if( z<0 || z>1)
		errorOut("mtLogisticR: opt.rFlag=1 and z should be [0,1]");
	else{
		for(i=0,j=0;i<m;i++){
			if(i >= ind[j+1])
				j++;
			b[i] = (p_flag[i] * m1[j]/(m1[j]+m2[j]))
				- ( (1^p_flag[i]) * m2[j]/(m1[j]+m2[j]));
		}
		
		// ATb = zeros(n,k);
		memset(ATb,0,sizeof(double)*n*k);

		// Compute ATb
		for(i=0;i<k;i++){
			normalizedmv_lowlevel('T',ind[i+1]-ind[i] ,n,&A[ind[i]],m,&b[ind[i]],tt,ind[i],0, &opt);
			dcopy(n,tt,1,&ATb[i*n],1);
		}

		if(q==1)
			q_bar = DBL_MAX;
		else if ( q> 1.0e6 )
			q_bar = 1.0;
		else
			q_bar = q/(q-1);
		lambda_max = 0.0;
		for(i=0; i<n ; i++)
			lambda_max = maxof2(lambda_max ,dnrm(k,&ATb[i],n,q_bar));
		lambda_max = lambda_max / m;
		lambda = z * lambda_max;
	}

	if( opt.init == 2){
		memset(x,0,sizeof(double)*n*k);
		memset(c,0,sizeof(double)*k);
	}
	else{
		if(hasGot(&opt,"x0"))
			x = opt.x0;
		else
			memset(x,0,sizeof(double)*n*k);
		if( hasGot(&opt,"mt_c0"))
			c = opt.mt_c0;
		else // c = log(m1./m2);
			for(i=0;i<k;i++)
				c[i] = log(m1[i]/m2[i]);
	}

	for(i=0; i < k; i++)
		normalizedmv_lowlevel('N',ind[i+1]-ind[i],n,&A[ind[i]],m,&x[i*n],&Ax[ind[i]],ind[i],0, &opt);

	if( opt.mFlag == 1 && opt.lFlag == 1 && opt.q==2){
		bFlag = 0;

		L = 1.0l/m;

		//xp=x; Axp=Ax; xxp=zeros(n,k);
		//cp=c; ccp=zeros(1,k);
		dcopy(n*k,x,1,xp,1);
		dcopy(m,Ax,1,Axp,1);
		memset(xxp,0,sizeof(double)*n*k);
		dcopy(k,c,1,cp,1);
		memset(ccp,0,sizeof(double)*k);

		//t=sqrt(sum(x.^2, 2)); tp=t;
		for(i=0; i<n ; i++)
			t[i] = dnrm2(k,&x[i],n);
		dcopy(n,t,1,tp,1);


		gamma=1;

		for(iterStep=0;iterStep<opt.maxIter;iterStep++){
			for(;;){
				if(iterStep!=0){
					alpha = (-gamma + sqrt(gamma*gamma + 4*L*gamma))/(2*L);
					beta = (gamma - gamma*alphap)/(alphap*gamma + alphap*L*alpha);
					// s = x + beta* xxp;
					dcopy(n*k,x,1,s,1);
					daxpy(n*k,beta,xxp,1,s,1);

					//sc = c + beta * ccp;
					dcopy(k,c,1,sc,1);
					daxpy(k,beta,ccp,1,sc,1);

					// s_t= t + beta * (t-tp);
					dcopy(n,t,1,s_t,1);
					dscal(n,1.0 + beta, s_t,1);
					daxpy(n,-beta,tp,1,s_t,1);

					// As = Ax + beta * (Ax - Axp);
					dcopy(m,Ax,1,As,1);
					dscal(m,1+beta,As,1);
					daxpy(m,-beta,Axp,1,As,1);
				}
				else{
					alpha = (-1 + sqrt(5))/2;
					beta = 0;
					//s=x; s_t=t; sc=c; As=Ax;
					dcopy(n*k,x,1,s,1);
					dcopy(n,t,1,s_t,1);
					dcopy(k,c,1,sc,1);
					dcopy(m,Ax,1,As,1);
				}

				dcopy(m,As,1,aa,1);
				//aa= - diag(y) * (A * s + sc)
				for(i=0; i < k; i++)
					daxpy(ind[i+1]-ind[i],sc[i],oneVector,1,&aa[ind[i]],1);
				dtbmv('U','N','N',m,0,y,1,aa,1);
				dscal(m,-1.0,aa,1);

				fun_s = 0;
				for(i=0; i<m;i++){
					bb[i] = maxof2(aa[i],0.0);
					fun_s += ( log(exp(-bb[i])+exp(aa[i]-bb[i])) + bb[i]);
				}
				fun_s = fun_s / m;

				//b = -y .* (1-prob);
				for(i=0;i<m;i++){
					prob[i] = 1.0l/(1.0l+exp(aa[i]));
					b[i] = (- y[i] * (1.0l - prob[i])) / m;
				}
				

				for(i=0;i<k; i++)
					gc[i] = ddot(ind[i+1]-ind[i],&b[ind[i]],1,oneVector,1);

				// Compute ATb
				for(i=0;i<k;i++){
					normalizedmv_lowlevel('T',ind[i+1]-ind[i] ,n,&A[ind[i]],m,&b[ind[i]],tt,ind[i],0, &opt);
					dcopy(n,tt,1,&g[i*n],1);
				}

				//u = s - g/L; 
				dcopy(n*k,s,1,u,1);
				daxpy(n*k,-1/L,g,1,u,1);

				//cnew = sc - gc/L;
				dcopy(k,sc,1,cnew,1);
				daxpy(k,-1/L,gc,1,cnew,1);

				//v= s_t - lambda / L;
				dcopy(n,s_t,1,v,1);
				daxpy(n,-lambda/L,oneVector,1,v,1);


				//eppMatrix(x,v,n,k,lambda/L,q);
				ep21R(xnew,tnew,u,v,n,k);

				//v = xnew-s;
				dcopy(n*k,xnew,1,v,1);
				daxpy(n*k,-1,s,1,v,1);

				// Compute Axnew
				for(i=0; i < k; i++)
					normalizedmv_lowlevel('N',ind[i+1]-ind[i],n,&A[ind[i]],m,&xnew[i*n],&Axnew[ind[i]],ind[i],0, &opt);

				dcopy(m,Axnew,1,aa,1);
				//aa= - diag(y) * (A * xnew + c)
				for(i=0; i < k; i++)
					daxpy(ind[i+1]-ind[i],cnew[i],oneVector,1,&aa[ind[i]],1);
				dtbmv('U','N','N',m,0,y,1,aa,1);
				dscal(m,-1.0,aa,1);

				fun_x = 0;
				for(i=0; i<m;i++){
					bb[i] = maxof2(aa[i],0.0);
					fun_x += ( log(exp(-bb[i])+exp(aa[i]-bb[i])) + bb[i]);
				}
				fun_x = fun_x / m;

				r_sum = (ddot(n*k,v,1,v,1) + ddot(k,cnew,1,cnew,1) 
					+ ddot(k,sc,1,sc,1) - 2*ddot(k,cnew,1,sc,1) -2*ddot(n,tnew,1,s_t,1) + ddot(n,tnew,1,tnew,1) + ddot(n,s_t,1,s_t,1))/2 ; 
				l_sum = fun_x - fun_s - ddot(n*k,v,1,g,1) 
					- ddot(k,cnew,1,gc,1) + ddot(k,sc,1,gc,1);

				if(r_sum <= 1e-20){
					bFlag=1;
					break;
				}

				if(l_sum <= r_sum * L)
					break;
				else
					L = maxof2( 2*L , l_sum/r_sum );
			}

			gamma = L * alpha * alpha;
			alphap = alpha;

			ValueL[iterStep] = L;

			tao = L * r_sum / l_sum;
			if(tao>=5)
				L = L * 0.8;

			//xp=x; x=xnew; xxp=x-xp;
			dcopy(n*k,x,1,xp,1);
			dcopy(n*k,xnew,1,x,1);
			dcopy(n*k,x,1,xxp,1);
			daxpy(n*k,-1,xp,1,xxp,1);

			//Axp=Ax; Ax=Axnew;
			dcopy(m,Ax,1,Axp,1);
			dcopy(m,Axnew,1,Ax,1);

			//tp=t; t=tnew;
			dcopy(n,t,1,tp,1);
			dcopy(n,tnew,1,t,1);

			//cp=c; c=cnew; ccp=c-cp;
			dcopy(k,c,1,cp,1);
			dcopy(k,cnew,1,c,1);
			dcopy(k,c,1,ccp,1);
			daxpy(k,-1,cp,1,ccp,1);


			fVal = fun_x + lambda*ddot(n,t,1,oneVector,1);

			fValp = funVal[iterStep] = fVal;

			if(iterStep!=0)
				fValp = funVal[iterStep-1];

			if(bFlag)
				break;
			norm_xp = dnrm2(n*k,xp,1);
			norm_xxp = dnrm2(n*k,xxp,1);
			if(terminationCondition(&opt,fVal,fValp,x,norm_xxp,norm_xp,iterStep))
				break;
		}
	}
	ret.errorCode = 0;
	ret.funVal = funVal;
	ret.totIter = iterStep;
	ret.ValueL = ValueL;
	ret.x = x;
	return ret;

}

